# Replication script for the Network Science article
# "Multimodal mechanisms of political discourse dynamics and the case of Germany’s nuclear energy phase-out"

# Import and setup data ---------

## Install and load packages -----
# rm(list=ls())
# install.packages(goldfish)
# remotes::install_github("snlab-ch/goldfish@v1.5.2")
library(goldfish)
library(tidyverse)

## Import main data -------------
### Import network data ---------
netdata <- readRDS(file = "akw_netdata.rds") %>% 
  rename(sender=speaker, receiver=claim, time=date, increment=weight) %>%
  select(sender, receiver, time, increment, everything()) %>%
  mutate(time = as.POSIXlt(time),
         govopp = coalesce(recode(party, 
                                  `31` = "government", `33` = "government",
                                  `32` = "opposition",`34` = "opposition",
                                  `35` = "opposition", 
                                  .default = NA_character_),
                           recode(sender, 
                                  `Bundesregierung` = "government", 
                                  .default = NA_character_))) %>%
  arrange(time) %>% as_tibble()
sum(netdata$increment==1)/nrow(netdata) # 78% (~3/4) positive claims
netdata

### Import and add power and political group attributes -----
actor_attributes <- read.csv2("actor_attributes.csv", stringsAsFactors = F)
netdata <- full_join(netdata, actor_attributes, by = c("sender" = "actor")) %>%
  mutate(category = coalesce(actor_category.x, actor_category.y),
         party = coalesce(party.x, party.y),
         org = coalesce(organization.x, organization.y),
         politician = ifelse(is.na(party), 0, 1),
         govt = coalesce(ifelse(party == 31 | party == 33, 1, 0), 0),
         opp = coalesce(ifelse(party == 32 | party == 34 | party == 35 | party == 36, 1, 0), 0),
         office = ifelse(category == 20, 1, 0)) %>%
  select(-c(actor_category.x, actor_category.y, party.x, party.y, organization.x, organization.y)) %>%
  ## Right-Left-Scale (rile) (2009-09)
  mutate(rile_pos = recode(party, `31` = 8.724, # CDU/CSU	8.724
                           `32` = -18.297, # SPD	-18.297	
                           `33` = 4.272, # FDP	4.272
                           `34` = -13.571, # Grüne	-13.571
                           `35` = -24.519, .default = NA_real_), # LINKE	-24.519
         coalition = recode(coalition, "pro" = 1, "contra" = 0, .default = NA_real_),
         eco_pos = recode(party, `31` = 5.833, # CDU/CSU	5.833
                          `32` = 5.254, # SPD	5.254	
                          `33` = 5.029, # FDP	5.029
                          `34` = 10.651, # Grüne	10.651
                          `35` = 5.529, .default = NA_real_) # LINKE	5.529
  ) %>% filter(!is.na(receiver))
netdata
## Party distances according to the Manifeso Project data
## Volkens, Andrea / Krause, Werner / Lehmann, Pola / Matthieß, Theres / 
## Merz, Nicolas / Regel, Sven / Weßels, Bernhard (2018): The Manifesto Data 
## Collection. Manifesto Project (MRG/CMP/MARPOR). Version 2018b. 
## Berlin: Wissenschaftszentrum Berlin für Sozialforschung (WZB). 
## https://doi.org/10.25522/manifesto.mpds.2018b

## Define the two node sets --------
### Speakers preexist --------------
speakers <- netdata %>% filter(!duplicated(sender)) %>%
  mutate(present = TRUE) %>%
  rename(label = sender, role = category) %>%
  select(label, present, org, party, role, power, politician, govt, opp, office,
         formal_political_actor, individual, coalition) %>% 
  defineNodes()
speakers

### Concepts raised at particular points -----
concepts <- data.frame(label=unique(netdata$receiver), present=FALSE,
                       stringsAsFactors = F) %>% defineNodes()
conceptstart <- netdata %>% filter(!duplicated(receiver)) %>%
  select(time, receiver) %>% rename(node=receiver) %>%
  mutate(time = as.POSIXct(time), replace=TRUE)
concepts <- linkEvents(concepts, changeEvents = conceptstart, attribute='present')
concepts

## Define main network variables --------
### Claim support ----------------------
posnet <- defineNetwork(nodes = speakers, nodes2 = concepts)
posevents <- netdata[,1:4] %>% filter(increment==1) %>%
  select(time, sender, receiver) %>%
  mutate(increment = 1)
posnet <- linkEvents(posnet, changeEvents = posevents, 
                     nodes = speakers, nodes2 = concepts)
posdep <- netdata %>%
  filter(increment==1 & formal_political_actor == 1) %>% 
  select(time, sender, receiver, increment)
posdep <- defineDependentEvents(posdep, 
                                speakers, concepts, 
                                defaultNetwork = posnet)

### Claim contestation ---------------
negnet <- defineNetwork(nodes = speakers, nodes2 = concepts)
negevents <- netdata[,1:4] %>% filter(increment==-1) %>%
  select(time, sender, receiver) %>%
  mutate(increment = 1)
negnet <- linkEvents(negnet, changeEvents = negevents, 
                     nodes = speakers, nodes2 = concepts)


## Construct periods of dependent events ----------
period1 <- netdata %>%
  filter(increment==1 & formal_political_actor == 1 &
           time >= "2011-03-11" & time < "2011-03-17") %>% 
  select(time, sender, receiver, increment)
period1 <- defineDependentEvents(period1, 
                                 speakers, concepts, 
                                 defaultNetwork = posnet)
period2 <- netdata %>%
  filter(increment==1 & formal_political_actor == 1 &
           time >= "2011-03-16" & time < "2011-04-09") %>% 
  select(time, sender, receiver, increment)
period2 <- defineDependentEvents(period2, 
                                 speakers, concepts, 
                                 defaultNetwork = posnet)
period3 <- netdata %>%
  filter(increment==1 & formal_political_actor == 1 &
           time >= "2011-04-09" & time <= "2011-05-30") %>% 
  select(time, sender, receiver, increment)
period3 <- defineDependentEvents(period3, 
                                 speakers, concepts, 
                                 defaultNetwork = posnet)
period4 <- netdata %>%
  filter(increment==1 & formal_political_actor == 1 &
           time >= "2011-05-31" & time <= "2011-07-01") %>% 
  select(time, sender, receiver, increment)
period4 <- defineDependentEvents(period4, 
                                 speakers, concepts, 
                                 defaultNetwork = posnet)
goldfishObjects()

# Estimation ---------
estPrefs <- list(returnIntervalLogL = T, initialDamping = 40, 
                 engine = "gather_compute", impute = T)

## Rate submodels -------------
### Period 1 -------------------
p1r <- estimate(period1 ~ 1 + outdeg(posnet) + 
                  ego(speakers$power) + ego(speakers$office) + ego(speakers$individual),
                model = "DyNAM", subModel = "rate",
                estimationInit = estPrefs, verbose = F)
tidy(p1r)
### Period 2 -------------------
p2r <- estimate(period2 ~ 1 + outdeg(posnet) + 
                  ego(speakers$power) + ego(speakers$office) + ego(speakers$individual),
                model = "DyNAM", subModel = "rate",
                estimationInit = estPrefs, verbose = F)
tidy(p2r)
### Period 3 -------------------
p3r <- estimate(period3 ~ 1 + outdeg(posnet) + 
                  ego(speakers$power) + ego(speakers$office) + ego(speakers$individual),
                model = "DyNAM", subModel = "rate",
                estimationInit = estPrefs, verbose = F)
tidy(p3r)
### Period 4 -------------------
p4r <- estimate(period4 ~ 1 + outdeg(posnet) +
                  ego(speakers$power) + ego(speakers$office) + ego(speakers$individual),
                model = "DyNAM", subModel = "rate",
                estimationInit = estPrefs, verbose = F)
tidy(p4r)

## Choice submodels ----------------
### Period 1 -------------------
p1c <- estimate(period1 ~ indeg(posnet) + indeg(negnet)
                + four(posnet) +
                  tertius(posnet, speakers$power, 
                          aggregateFun = function(x) max(x, na.rm = TRUE)) +
                  tertius(posnet, speakers$party, 
                          aggregateFun = function(x) 
                            vegan::diversity(table(as.character(na.omit(x))))) +
                  tertiusDiff(posnet, speakers$govt, 
                              aggregateFun = function(x) mean(x, na.rm = TRUE)),
                model = "DyNAM", subModel = "choice",
                estimationInit = estPrefs, verbose = F)
tidy(p1c)

### Period 2 -------------------
p2c <- estimate(period2 ~ indeg(posnet) + indeg(negnet)
                + four(posnet) +
                  tertius(posnet, speakers$power, 
                          aggregateFun = function(x) max(x, na.rm = TRUE)) +
                  tertius(posnet, speakers$party, 
                          aggregateFun = function(x) 
                            vegan::diversity(table(as.character(na.omit(x))))) +
                  tertiusDiff(posnet, speakers$govt, 
                              aggregateFun = function(x) mean(x, na.rm = TRUE)),
                model = "DyNAM", subModel = "choice",
                estimationInit = estPrefs, verbose = F)
tidy(p2c)

### Period 3 -------------------
p3c <- estimate(period3 ~ indeg(posnet) + indeg(negnet)
                + four(posnet) +
                  tertius(posnet, speakers$power, 
                          aggregateFun = function(x) max(x, na.rm = TRUE)) +
                  tertius(posnet, speakers$party, 
                          aggregateFun = function(x) 
                            vegan::diversity(table(as.character(na.omit(x))))) +
                  tertiusDiff(posnet, speakers$govt, 
                              aggregateFun = function(x) mean(x, na.rm = TRUE)),
                model = "DyNAM", subModel = "choice",
                estimationInit = estPrefs, verbose = F)
tidy(p3c)

### Period 4 -------------------
p4c <- estimate(period4 ~ indeg(posnet) + indeg(negnet)
                + four(posnet) +
                  tertius(posnet, speakers$power, 
                          aggregateFun = function(x) max(x, na.rm = TRUE)) +
                  tertius(posnet, speakers$party, 
                          aggregateFun = function(x) 
                            vegan::diversity(table(as.character(na.omit(x))))) +
                  tertiusDiff(posnet, speakers$govt, 
                              aggregateFun = function(x) mean(x, na.rm = TRUE)),
                model = "DyNAM", subModel = "choice",
                estimationInit = estPrefs, verbose = F)
tidy(p4c)

# Diagnostics -------------

## Check outliers ----------
examineOutliers(p1r)
examineOutliers(p2r)
examineOutliers(p3r)
examineOutliers(p4r)
examineOutliers(p1c)
examineOutliers(p2c)
examineOutliers(p3c)
examineOutliers(p4c)

## Check changepoints ----------
examineChangepoints(p1r)
examineChangepoints(p2r)
examineChangepoints(p3r)
examineChangepoints(p4r)
examineChangepoints(p1c)
examineChangepoints(p2c)
examineChangepoints(p3c)
examineChangepoints(p4c)

## Check whether a single period has changepoints ---------
p0r <- estimate(posdep ~ 1 + outdeg(posnet) + 
                  ego(speakers$individual) + ego(speakers$office) + 
                  ego(speakers$power),
                model = "DyNAM", subModel = "rate",
                estimationInit = estPrefs, verbose = F)
tidy(p0r)
examineChangepoints(p0r)
examineOutliers(p0r)
p0c <- estimate(posdep ~ indeg(posnet) + indeg(negnet) + 
                  four(posnet) +
                  tertius(posnet, speakers$power, 
                          aggregateFun = function(x) max(x, na.rm = TRUE)) +
                  tertius(posnet, speakers$party, 
                          aggregateFun = function(x) 
                            vegan::diversity(table(as.character(na.omit(x))))) +
                  tertiusDiff(posnet, speakers$govt, 
                              aggregateFun = function(x) sum(x, na.rm = TRUE)),
                model = "DyNAM", subModel = "choice",
                estimationInit = estPrefs, verbose = F)
tidy(p0c)
examineChangepoints(p0c)
examineOutliers(p0c)
